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The  preoccupation  of  geodesists  with  direct  and  inverse  position  computa¬ 
tion  on  the  terrestrial  ellipsoid  has  a  long  and  distinguished  history. 
Gauss,  Bessel,  Helmert,  Puissant,  Rainsford,  McCaw,  and  Sodano  are  all 
prominent  names  associated  with  formulas  and  algorithms  developed  for  the 
solution  of  what  in  the  German  technical  literature  used  to  be  called  "die 
geodaetische  Hauptaufgabe"  (the  principal  geodetic  problem).  With  the  advent 
of  electronic  computers,  this  work  of  giants  was  given  a  capstone  by  T. 
Vincenty  with  his  optimal  adaptation  for  automatic  computation  of  the  globally 
accurate  Bessel-Helmert-Rainsford  iterative  algorithms  (Vincenty  1975,  1976). 

The  "direct"  problem  can  be  posed  as  follows:  Given  the  position 
(latitude  and  longitude)  of  a  point  on  the  reference  ellipsoid  (the 
"standpoint"),  as  well  as  the  orientation  (forward  azimuth)  and  length  of  a 
geodesic  line  emanating  from  it,  compute  the  position  (latitude  and  longitude) 
of  the  terminal  point  of  that  geodesic  line  (the  "forepoint")  and  its  back 
azimuth.  The  "inverse"  problem,  as  can  be  expected,  is  the  converse  of  the 
direct  problem:  Given  the  coordinates  of  two  points  on  the  reference 
ellipsoid,  compute  the  length  of  the  geodesic  line  joining  them,  as  well  as 
the  forward  and  back  azimuths  at  the  respective  endpoints,  which  in  this  case 
are  arbitrarily  taken  to  be  the  standpoint  and  the  forepoint. 

Vincenty' s  direct  and  inverse  position  computation  algorithms  are 
efficient  and  accurate  to  a  fraction  of  a  millimeter  for  short  and  long  geode¬ 
sics  alike,  ranging  in  length  from  a  few  centimeters  to  just  under  half-way 
around  the  world.  As  such,  they  are  yardsticks  against  which  the  performance 
of  other  direct  and  inverse  position  computation  algorithms  are  to  be 
measured.  However,  they  are  iterative,  which  is  to  say  that  the  number  of 
steps  in  a  solution  can  vary  depending  on  the  geometry  of  the  individual 
problem. 

Since  convergence  is  very  fast  (two  or  three  iterations  is  the  norm),  the 
iterative  nature  of  Vincenty 'a  algorithms  is  hardly  a  consideration  in  non¬ 
realtime  applications  run  on  present-day  powerful  minicomputers  and 
mainframes;  in  faot,  it  has  been  shown  that  Vincenty 's  algorithms  execute 
faster  than  either  Sodano 's  or  Andoyer-Lambert's  noniterative  long-line 
counterparts.  There  are,  however,  applications  for  whioh  one  would 
intuitively  prefer  the  noniterative  solutions  of  the  dlreet  and  inverse 
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position  computation  problems,  solutions  which  would  not  have  the  complexity 
of  the  existing  long-line  noniterative  algorithms,  result  in  even  more  compact 
code,  and  execute  faster  than  Vincenty's  algorithms,  and  still  deliver  the 
desired  accuracy. 

One  such  application  is  the  computation  of  geodetic  survey  work  on  the 
ellipsoid  (as  opposed  to  plane  coordinates)  implemented  on  a  portable  micro¬ 
computer  (e.g.,  the  surveyor's  field  computer) .  Here  one  is  typically  faced 
with  very  limited  memory  and  the  need  for  compact  code,  with  speed  of  execu¬ 
tion  being  an  important  but  secondary  consideration.  Sub-millimeter 
computational  accuracy  is  required  in  this  application;  however,  the  line 
length  is  limited  by  intervisibility  and  seldom  exceeds  50  km. 

Another  such  application  occurs  in  hydrographic  surveying;  i.e.,  the 
realtime  computation  of  the  position  of  a  survey  vessel  with  respect  to  shore 
control  stations,  when  one  wishes  to  work  with  the  reference  ellipsoid  rather 
than  with  a  map  projection.  Here  the  maximum  line  length  will  vary  from  50  km 
for  line-of-sight  positioning  systems,  to  300  km  for  medium-range  systems  such 
as  Raydist  or  Argo,  to  1500  km  for  a  long-range  system  such  as  LORAN-C.  On 
the  other  hand,  the  computational  accuracy  requirement  can  be  proportionately 
relaxed  by  two,  three,  and  four  orders  of  magnitude  (compared  with  the 
geodetic  surveying  case),  depending  on  the  scale  of  the  survey  and  positioning 
method  used;  e.g.,  0.1  m  for  short-range  control,  large-scale  surveys  (for 
harbor  approach  charts),  1  m  for  medium-range  control,  medium-scale  surveys 
(for  coastal  sailing  charts),  and  10  m  for  long-range  control,  small-scale 
surveys  (for  general  sailing  charts).  Since  in  a  realtime  application  the 
respective  algorithms  must  execute  within  an  assigned  time  slot  measured  in 
milliseconds,  speed  of  execution  is  the  primary  consideration  in  this 
instance. 

Recently,  B.  R.  Bowring  of  Surrey,  England,  developed  and  published  very 
elegant  noniterative  algorithms  for  the  direct  and  Inverse  position  computa¬ 
tion  over  "short"  geodesic  lines  up  to  150  km  (Bowring  1981).  These  "quasi- 
spherioal"  formulas  are  remarkably  concise  and  accurate  within  their  intended 
range  of  application;  they  very  likely  represent  the  last  word  in  streamlining 
the  solution  of  the  "principal  geodetic  problem."  Bowring's  algorithms  lend 
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themselves  admirably  to  the  first  application  outlined  above,  i.e. ,  computa¬ 
tion  of  geodetic  survey  work.  They  were  successfully  used  by  this  writer  as 
the  basis  of  a  powerful  and  efficient  geodetic  package  of  programs  implemented 
on  the  Hewlett-Packard  HP-9815A  desktop  computer  (Taylor  1981). 

The  purpose  of  this  paper  is  to  present  the  results  of  an  investigation 
as  to  the  extent  to  which  Bowring's  algorithms  are  sufficiently  accurate  to 
support  the  second  application  mentioned  above,  i.e.,  the  realtime  positioning 
of  a  surface  vessel  for  hydrographic  surveying  or  precise  navigation  purposes. 
This  investigation  evaluated  the  total  position  error  produced  by  Bowring's 
algorithms  over  a  large  number  of  geodesic  lines  emanating  from  standpoints 
located  at  seven  representative  latitudes  (0,  15,  30,  45,  60,  75,  and  89 
degrees),  in  nine  representative  azimuths  (0,  30,  45,  60,  90,  120,  135,  150, 
and  180  degrees),  and  of  lengths  ranging  from  50  to  4000  km  (preliminary  com¬ 
putations  indicated  this  distance  to  be  the  usable  limit).  In  all,  4284  cases 
were  computed. 

As  a  first  step  in  every  case,  the  coordinates  of  each  forepoint  were 
computed  using  the  precise  Vineenty's  direct  algorithm.  This  step  was  then 
repeated  using  Bowring's  direct  algorithm,  and  the  distance  separating  the  two 
sets  of  coordinates  was  taken  as  the  total  position  error  of  Bowring's  direct 
algorithm.  Next,  Bowring's  inverse  algorithm  was  used  to  recover  the  length 
and  forward/back  azimuth  of  the  geodesic  line  between  the  given  standpoint  and 
computed  precise  forepoint  coordinates.  The  resulting  length  and  forward 
azimuth  were  then  used  as  arguments  in  Vineenty’s  direct  algorithm  to  compute 
another  set  of  forepoint  coordinates,  and  the  distance  separating  the  two  sets 
of  coordinates  was  taken  as  the  total  position  error  of  Bowring's  inverse 
algorithm. 

The  total  position  errors,  obtained  in  meters,  were  shown  as  proportional 
errors  relative  to  the  length  of  the  geodesic  line,  in  parts  per  million 
(ppm).  Inspection  of  the  tabulations  confirmed  that  Bowring's  direct  and 
inverse  algorithms  are  well  balanced  with  respect  to  accuracy,  as  the  corres¬ 
ponding  errors  in  any  given  case  were  always  very  nearly  equal.  Por  eaoh  of 
the  63  geodesic  lines  computed  at  distanoe  Increments  from  50  to  4000  km  (9 
radial  lines  from  eaoh  of  7  standpoint  latitudes),  the  tabulations  were 
searched  to  determine  the  distances  at  whioh  total  position  errors 


exceed  the  thresholds  of  0.0001  m,  0.001  m,  0.01  m,  0.1m,  1  m,  10  m,  and  100 
m;  and  in  terms  of  relative  error,  the  thresholds  of  0.1  ppm  (1:10,000,000), 
0.2  ppm  (1:5,000,000),  1  ppm  (1:1,000,000),  2  ppm  (1:500,000),  10  ppm 

(1:100,000),  and  20  ppm  (1:50,000). 

As  could  be  expected  from  the  nature  of  the  problem,  the  worst  perfor¬ 
mance  for  radial  lines  emanating  from  each  standpoint  was  along  the  meridian, 
with  progressively  better  performance  along  geodesics  in  azimuths  away  from 
the  meridian.  For  each  of  the  seven  standpoint  latitudes,  this  worst-case 
performance  was  taken  as  the  upper  bound  of  the  total  position  error  to  be 
expected  of  Bowring's  direct  and  inverse  algorithms  over  any  geodesic  line 
having  an  endpoint  at  that  latitude.  The  resulting  information  is  portrayed 
graphically  in  Figures  1  and  2. 

By  inspection  of  the  log-linear  graph  of  Figure  1,  it  is  clear  that  even 
in  the  worst  possible  case  (geodesic  line  on  or  near  the  meridian  originating 
at  or  near  the  latitude  of  45  degrees),  Bowring's  algorithms  meet  the  compu¬ 
tational  accuracy  requirements  of  both  geodetic  survey  work  and  of  surface 
vessel  position  fixing  consistent  with  the  accuracy  of  shore-based  positioning 
systems  likely  to  be  used  for  hydrographic  surveying  and  precise  navigation 
purposes.  The  total  position  error  produced  by  either  the  direct  or  the 
inverse  algorithm  is  guaranteed  to  be  less  than  0.001  m  up  to  100  km,  less 
than  0.1  m  up  to  500  km,  and  less  than  10  m  up  to  1500  km.  One  notes  that  the 
error  curves  are  symmetrical  about  the  latitude  of  45  degrees,  and  that 
progressively  better  accuracy  performance  is  obtained  along  geodesic  lines 
originating  in  both  lower  and  higher  latitudes,  as  well  as  along  geodesics  in 
azimuths  away  from  the  meridian. 

It  is  also  interesting  to  note  that  on  the  log-log  graph  of  Figure  2,  the 
relative  error  in  parts  per  million  as  a  function  of  line  length  is  linear. 
This  quite  unexpected  result  clearly  suggests  an  empirical  formula  for  the 
global  upper  bound  of  the  total  position  error  produced  by  Bowring's 
algorithms.  By  considering  the  worst-case  performance,  the  following 

empirical  relationship  (Equation  (1))  for  the  maximum  relative  error  in  parts 
per  million  (Mppn)  as  a  funotlon  of  geodesic  line  length  in  kilometers  (Dfem) 
can  be  derived: 


Mppa«  7.17x10-9  Dg*86 


(1) 
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An  expression  for  the  absolute  maximum  error  can  be  derived  by  multiplying 
Equation  (1)  by  distance.  Taking  into  account  the  conversion  of  distance 
units  to  meters,  the  following  equation  results: 

Mb  =  7.17x10-12  0^86  (2) 

The  error  curves  in  Figure  1  depict  somewhat  more  conservative  error  estimates 
than  the  values  given  by  Equation  (2).  This  is  due  to  upward  rounding  of  the 
total  position  errors  on  the  computer  printout  from  which  data  shown  in  Figure 
1  were  compiled. 

Listings  of  a  FORTRAN  implementation  of  Bowring's  direct  and  inverse 
algorithms  are  given  in  Figures  3  and  4,  and  those  of  Vincenty's  direct  and 
inverse  algorithms  in  Figures  5  and  6. 
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